Hypotheses

H1. Forest structure is negatively correlated with increasing elevation. H2. Forest structure is negatively correlated with increasing slope angle. H3. Forest structure is negatively correlated with north and east aspects. H4. Forest structure is controlled/determined by topographic position index/ landform classifications.
(forest structural homogeneity is positively correlated with increasing spatial resolution of the topographic position index) H5. Forest structure is negatively correlated with increasing Topographic Roughness. H6. Forest structure will be positively correlated with land cover and forest type.

Site Description

The location chosen for this project is Bighorn National Forest in Wyoming, United States. Bighorn National Forest is approximately 128 meters long and 48 meters wide, encompassing over 404685 hectares of land (USDA Forest Service 2021). The elevation of Bighorn National Forest ranges from 1524 meters, to its highest elevation of 4013 meters at Cloud Peak (USDA Forest Service 2021). Approximately 66% of the land cover within Bighorn National Forest is forested and made up of Douglas-fir, ponderosa pine, lodgepole pine, and spruce (Witt 2008). The sites chosen for this study are located within the National Elevation Dataset N35 W108, as shown with the black outline in Figure X, which is approximately 4828 square kilometers. Bighorn National Forest

Site Selection

From August 12-21, 2021, seventy-four plots at 10x10 meters were selected across the study area of Bighorn National Forest. The sites chosen were representative of the adjacent area, with their suitability primarily derived by their elevation to ensure they were within the montane biogeographic zone (between 1700 and 3000 meters). [Elevations lower than 1700 meters, and greater than 3000 meters were derived from a 1/3 arc-second digital elevation model (DEM) were contained in polygons and excluded from site selection consideration. The secondary consideration for site selection was given to the accessibility of the site. The majority of the sites chosen were on average 100 meters from a secondary road (dirt or gravel and not maintained), with few sites being greater than 100 meters from a primary road (paved and maintained). Topographic and geographic features were also included in the analysis of site accessibility, with steep inclines or declines, cliffs, major water body crossings and unstable or hazardous features generally excluded from consideration. Sites that had obvious signs of logging activity or forest fires were also excluded from consideration. The last criteria for potential sites were that they be analogous to their surrounding area, and contained at least one mature tree species (Diameter Breast Height (DBH) must be 10 centimeters or greater). After the previous criteria were met, final site selection was derived stochastically by walking a number of paces in one direction determined by a random number generator. The location produced from the previous step was used as the Northwest coordinate of the plot. Bighorn National Forest

Field Data:

For each site, a 10m x 10m plot was demarcated. The plots were aligned north-south, east-west, with markers placed at the northwest, northeast, southeast and southwest corners. The preliminary elevation and photos of the plot were taken at the northwest corner (Fig. 3). At each corner, as well as at the center of each plot, a CI-110 Plant Canopy Imager recorded the coordinates, as and a wide-angle image of the plot canopy. Additionally, the Plant Canopy Imager generated canopy data for the plot including the Gap Fraction, Photosynthetically Active Radiation (PAR) Leaf Area Index (LAI), PAR Average, Sunflecks, Canopy Density, and Leaf Angle. Geographic and topographic data about each plot was noted as well as the presence or absence of rocks, boulders, downed woody debris (DWD) or standing dead trees (snags), water and ground junipers. The vegetative species composition of the plot was identified and for trees with a DBH of 10cm or greater, DBH, richness and occurrence was documented.

Plot Photos

Libraries

library(moments)
library(gvlma)
library(raster)
library(lmtest)
library(car)
library(leaps)
library(tidyverse)
library(leaflet)
library(ggplot2)
library(ggthemes)
library(tidyr)
library(broom)
library(sp)
library(sf)
library(vegan)
library(dplyr)
library(spData)
library(stats)
library(evaluate)
library(spatialEco)
library(knitr)
library(insol)
library(tmap)
library(geodiv)
library(GmAMisc)
library(viridis)
library(rgdal)
library(lidR)
library(lattice)
library(rasterVis)
library(RColorBrewer)
library(ggcorrplot)
library(cowplot)
library(googleway)
library(ggrepel)
library(ggspatial)
library(rnaturalearth)
library(rnaturalearthdata)
library(citation)
library(MASS)
library(lme4)

knitr::opts_chunk$set(cache=TRUE, warning=FALSE, message=FALSE)

Read in Field Data

CanopyData <- read.csv("proj_data/CanopyData.csv")
PlotData <- read.csv("proj_data/PlotData.csv")
TreeData <- read.csv("proj_data/TreeData.csv")
TreeHeightsandAge <-read.csv("proj_data/Plot_tree_stats.csv")
DBHsAndDensity <-read.csv("proj_data/DBH_Density.csv")

Bighorn National Forest Shape File

nf <- read_sf("proj_data/S_USA.AdministrativeForest/S_USA.AdministrativeForest.shp") %>%
  st_transform(.,crs=32613)

bighorn <- dplyr::filter(nf,FORESTNAME=="Bighorn National Forest")

Load DEM Data

DEM <- raster("proj_data/DEM.tif")

Elevation of Bighorn National Forest

Elevation of Bighorn National Forest.

Load in Wyoming Landuse Landcover Data

#Read in the LULC Raster 
lulc=raster("proj_data/lulc.tif")

# LULC Names 
lulc_names =read_csv("proj_data/LULC.csv") %>%
  dplyr::select(class=ECOLSYS_LU...2, ID=Value)

# Convert to a Raster
lulc=as.factor(lulc)

# Update the RAT with a Left Join
levels(lulc)=left_join(levels(lulc)[[1]],lulc_names)
#table(values(lulc))

Plot of Land Use Land Cover

lulc_colors=data.frame(class=levels(lulc)[[1]]$class)

lulc_colors$col=c("#EBA998","#FFFB00", "#4080BF", "#989F60", "#AF5063", "#3200FF", "#A289A8", "#BE20DF", "#971452", "#6250AF" ,  "#8EB880", "#624074", "#A0B232", "#B232A0", "#FF00E8","#F1197C", "#7C2B4B", "#7C2B72", "#A0D092", "#D095DE", "#70808F", "#10D4EF", "#AF50A7", "#EF103D","#4F30CF", "#FFC300", "#BCFF00", "#00BBFF", "#2080DF", "#DF8E20", "#9A8A86", "#CF3050", "#708F75", "#8F708D", "#95BCDE",  "#8F8F70", "#FFCD00", "#EB98E1", "#543E6C", "#DFBA20", "#FF9300", "#00FF43", "#B25932", "#B2327F", "#543B1B", "#1B543B", "#1B3D54", "#805C73", "#B62B85", "#F119B2", "#F18E19", "#591AAF", "#E0FF00", "#E0FF90")

lulcplot<-rasterVis::gplot(lulc)+
  geom_raster(aes(fill=as.factor(value)))+
  scale_fill_manual(labels=lulc_colors$class,
                    values=lulc_colors$col,
                    name="Landcover Type")+
  coord_equal()+
  theme(legend.position = "bottom", legend.key.size=unit(.5, "line"), legend.text=element_text(size=8))+
  guides(fill=guide_legend(ncol=3,byrow=TRUE))

lulcplot

Rename, Summarize, and Group Data

canopy <- CanopyData %>%
  group_by(Plot_Number) %>%
  dplyr::select(-Sunflecks) %>% 
  summarise_if(is.numeric, list(mean = function(x) mean(x, na.rm = TRUE), 
                                      sd = function(x) sd(x, na.rm = TRUE)))
plots <- PlotData %>%
 dplyr::group_by(Plot_Number) %>%
 dplyr::summarize(lat=mean(Lat), lon=mean(Lon)) 

trees<- TreeData %>%
  summarize(Plot_Number, Tree_ID, Species_Scientific_Name, DBH, Aprox_Tree_Age) %>%
   dplyr:: group_by(Plot_Number)

Study Site Plots

#GGPlot of Study Site
studysiteplots<-ggplot(plots)+
  geom_point(mapping=aes(x=lon, y=lat, color=Plot_Number)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Longitude", y ="Latitude",
       title ="Study Site Plots", color ="Plot Number") + 
  theme(plot.title = element_text(hjust = 0.5))
studysiteplots

#Isolate standard lidar study plots 
Lidarplots <-plots[-c(7,12,13, 23,24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 43,
                              51,57, 58,59, 65,66, 69),]

#GGPlot of Lidar Plots 
lidarplotsgg<-ggplot(Lidarplots)+
  geom_point(mapping=aes(x=lon, y=lat, color=Plot_Number))+
  scale_color_viridis_c(option="turbo")+
       labs(x = "Longitude", y ="Latitude",
       title ="Lidar Plots", color ="Plot Number") + 
  theme(plot.title = element_text(hjust = 0.5))
lidarplotsgg

Table of Species by Plot

tree1<-trees %>% 
  dplyr::select(-Tree_ID, -DBH) %>% 
  dplyr::group_by(Plot_Number,Species_Scientific_Name) %>% 
  dplyr::summarise(count=n()) %>%
  tidyr::spread(Species_Scientific_Name,count) %>% 
  dplyr::mutate_all(function(x) ifelse(is.na(x),0,x))

treecounts <-table(trees$Species_Scientific_Name)
view(tree1)
view(treecounts)

tree_coords <- merge(trees, plots)

Trees

This graph displays the location of species across plots.

spatial_tc<-ggplot(data=tree_coords, mapping=aes(x=lon, y=lat, color=Species_Scientific_Name))+
  geom_point()+  
  scale_color_viridis_d(option="turbo")+
  labs(x = "Latitude", y = "Longitude", title ="Spatial Occurrence of Tree Species", color ="Species Scientific Name") +
  theme(plot.title = element_text(hjust = 0.5))

spatial_tc

Species Occurence Throughout Plots

This bar graph displays the total observations of each species across all plots. The Pinus contorta or Lodgepole Pine was observed the most frequently with more than 350 sightings. The Pinus flexilis, Limber Pine, and the Pinus ponderosa, the Ponderosa Pine each only had one sighting across the plots.

tree_counts_bargraph<-barplot(treecounts,
        main= "Species Occurence Throughout Plots", 
        horiz=TRUE, 
        xlab="Occurence", 
        ylab= "Species",
        col = c("#DD8317", "#3A9B44", "#47ACC2", "#EACF4F", "#EA5F4f", "#DD1794","#FFC300"), 
        names.arg=c("Juniperus occidentalis", "Picea engelmannii","Pinus contorta", "Pinus ponderosa", "Pinus flexilis", "Populus tremuloides", "Pseudotsuga menziesii"),
        cex.names=0.2, legend=TRUE)

Calculate TPI, Slope angle, Slope aspect and Elevation

tpi10m <- raster("proj_data/tpi10m.tif")
tpi100m <- raster("proj_data/tpi100m.tif")
tpi1000m <- raster("proj_data/tpi1000m.tif")


TRI <- raster("proj_data/TRI.tif")
slope <- raster("proj_data/slope.tif")

Slope Angle of Bighorn National Forest.

aspect <- raster("proj_data/aspect.tif")

Slope Aspect of Bighorn National Forest.

eastwest <- raster("proj_data/eastwest.tif")
northsouth <- raster("proj_data/northsouth.tif")

Create Landforms at 10 Meter Resolution

SD10m <- sd(tpi10m[],na.rm=T) 
SD10m
## [1] 0.4458888
landform_sd10m <- reclassify(tpi10m, matrix(c(-Inf, -SD10m, 1,
                                              -SD10m, -SD10m/2, 2,
                                              -SD10m/2, 0, 3,
                                              0, SD10m/2, 4,
                                              SD10m/2, SD10m, 5,
                                              SD10m, Inf, 6),
                                            ncol = 3, byrow = T),
                             right = T)

# Turn it into categorical raster
landform_sd10m <- as.factor(landform_sd10m) 

rat_sd10m <- levels(landform_sd10m)[[1]]

rat_sd10m[["landform_sd10m"]] <- c('Valley', 'Lower Slope', 
                                   'Flat Area','Middle Slope', 
                                   'Upper Slope', 'Ridge')

levels(landform_sd10m) <- rat_sd10m


# Plot the classification
tpi10mplot <-rasterVis::levelplot(landform_sd10m, col.regions = rev(brewer.pal(6,'RdYlBu')),
                       labels = rat_sd10m$landcover,
                       main = "10M TPI Landform Classification",
                       colorkey=list(labels=list(at=1:6, labels=rat_sd10m[["landform_sd10m"]])))

tpi10mplot

histtpi10m <-hist(tpi10m[])

Create Landforms at 100 Meter Resolution

SD100m <- sd(tpi100m[],na.rm=T)
SD100m 
## [1] 2.491715
# Make Landform Classifications 
#Morphologic class De Reu et al. 2013;  Weiss (2001)
landform_sd100m <- reclassify(tpi100m, matrix(c(-Inf, -SD100m, 1,
                                                -SD100m, -SD100m/2, 2,
                                                -SD100m/2, 0, 3,
                                                0, SD100m/2, 4,
                                                SD100m/2, SD100m, 5,
                                                SD100m, Inf, 6),
                                              ncol = 3, byrow = T),
                              right = T)

# Turn it into categorical raster
landform_sd100m <- as.factor(landform_sd100m) 

rat_sd100m <- levels(landform_sd100m)[[1]]

rat_sd100m[["landform_sd100m"]] <- c('Valley', 'Lower Slope', 
                                     'Flat Area','Middle Slope', 
                                     'Upper Slope', 'Ridge')

levels(landform_sd100m) <- rat_sd100m
#writeRaster(landform_sd100m, file="proj_data/landform_sd100m.grd", overwrite=TRUE)

# Plot the classification
tpi100mplot<- rasterVis::levelplot(landform_sd100m, col.regions = rev(brewer.pal(6,'RdYlBu')),
                        labels = rat_sd100m$landcover,
                        main = "100m TPI Landform Classification",
                        colorkey=list(labels=list(at=1:6, labels=rat_sd100m[["landform_sd100m"]])))


tpi100mplot

histtpi100m <- histtpi100<-hist(tpi100m[])

Create Landforms at 1000 Meter Resolution

SD1000m <- sd(tpi1000m[],na.rm=T)
SD1000m 
## [1] 29.7036
# Make landform classes
#Morphologic class De Reu et al. 2013;  Weiss (2001)
landform_sd1000m <- reclassify(tpi1000m, matrix(c(-Inf, -SD1000m, 1,
                                                  -SD1000m, -SD1000m/2, 2,
                                                  -SD1000m/2, 0, 3,
                                                  0, SD1000m/2, 4,
                                                  SD1000m/2, SD1000m, 5,
                                                  SD1000m, Inf, 6),
                                                ncol = 3, byrow = T),
                               right = T)

# Turn it into categorical raster
landform_sd1000m <- as.factor(landform_sd1000m) 

rat_sd1000m <- levels(landform_sd1000m)[[1]]

rat_sd1000m[["landform_sd1000m"]] <- c('Valley', 'Lower Slope', 
                                       'Flat Area','Middle Slope', 
                                       'Upper Slope', 'Ridge')

levels(landform_sd1000m) <- rat_sd1000m 
#writeRaster(landform_sd1000m, file="proj_data/landform_sd1000m.grd", overwrite= TRUE) #tif wont work for as.factor 

# Plot the classification
tpi1000mplot<- levelplot(landform_sd1000m, col.regions = rev(brewer.pal(6,'RdYlBu')),
                         labels = rat_sd1000m$landcover,
                         main = "1000m TPI Landform Classification",
                         colorkey=list(labels=list(at=1:6, labels=rat_sd1000m[["landform_sd1000m"]])))

tpi1000mplot

histtpi1000m<-hist(tpi1000m[])

Lidar Images of Plots

#### Read in Plots
plotsf<- plots %>%
  st_as_sf(coords = c(x="lon", y="lat")) %>%
  st_set_crs(.,value=4326) 

# Obtain Plot Coordinates in Lidar Projection
lasfile="proj_data/LAS/plot_52_USGS_LPC_WY_Sheridan_2020_D20_13TCK1344.laz"

lproj=readLASheader(lasfile) %>% 
  crs()

#10m polys 
plot_polys <- plotsf %>% 
  st_transform(.,crs="+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m no_defs") %>% 
  st_buffer(dist = 50,endCapStyle="SQUARE") %>% 
  st_transform(.,crs=st_crs(4326)) 

# Plot Specific Change 
plotgeo <-plot_polys[52,] %>% 
  st_transform(.,crs=lproj) 


# Read in New LAS File
lasPlots <-readLAS(lasfile,
                   filter = paste("-keep_xy ",paste(st_bbox(plotgeo),
                                                    collapse=" "))) %>% 
  clip_roi(geometry= plotgeo) 

plot(lasPlots)

Plot 2.1, Lidar Image of Plot 2 in Bighorn National Forest. Plot 2.2, Lidar Image of Plot 2 in Bighorn National Forest. Plot 6, Lidar Image of Plot 6 in Bighorn National Forest. Plot 8.1, Lidar Image of Plot 8 in Bighorn National Forest. Plot 8.2, Lidar Image of Plot 8 in Bighorn National Forest. Plot 11.1, Lidar Image of Plot 11 in Bighorn National Forest. Plot 11.2, Lidar Image of Plot 11 in Bighorn National Forest. Plot 14.1, Lidar Image of Plot 14 in Bighorn National Forest. Plot 14.2, Lidar Image of Plot 14 in Bighorn National Forest. Plot 15, Lidar Image of Plot 15 in Bighorn National Forest. Plot 17, Lidar Image of Plot 17 in Bighorn National Forest. Plot 18, Lidar Image of Plot 18 in Bighorn National Forest. Plot 20, Lidar Image of Plot 20 in Bighorn National Forest. Plot 42.1, Lidar Image of Plot 42 in Bighorn National Forest. Plot 42.2, Lidar Image of Plot 42 in Bighorn National Forest. Plot 46, Lidar Image of Plot 46 in Bighorn National Forest. Plot 47.1, Lidar Image of Plot 47 in Bighorn National Forest. Plot 47.2, Lidar Image of Plot 47 at 30m Resolution in Bighorn National Forest. Plot 48.1, Lidar Image of Plot 48 at 30m Resolution in Bighorn National Forest. Plot 48.2, Lidar Image of Plot 48 at 30m Resolution in Bighorn National Forest. Plot 50, Lidar Image of Plot 50 in Bighorn National Forest. Plot 52.1, Lidar Image of Plot 52 at 100m Resolution in Bighorn National Forest. Plot 52.2, Lidar Image of Plot 52 at 30m Resolution in Bighorn National Forest.

Calculating Tree Heights

#dtm, dsm, chm 
dsm2 <- lidR::grid_canopy(lasPlots, res=2, p2r()) 
dtm2 <- lidR::grid_terrain(lasPlots, res= 2, algorithm=tin()) 


#Crop! 
chm2<-dsm2-dtm2 

chmcrop <- crop(chm2, plotgeo)  
plot(chmcrop)

#TRUE CROPPED MEAN 
chmmean<- cellStats(chmcrop, stat=mean) 
view(chmmean)  

#TRUE CROPPED SD
chmstd<-  cellStats(chmcrop, stat=sd) 
view(chmstd) 

#Max Value
chmmax <- cellStats(chmcrop, stat=max)
view(chmmax)

#Min Value 
chmmin <- cellStats(chmcrop, stat=min)
view(chmmin)

#Plot and Hist
plot(chmcrop, col = height.colors(25))

hist(chmcrop)

Calculate the Diversity

diversity=data.frame(Plot_Number=tree1$Plot_Number,div=diversity(tree1))

Species Diversity

This bar graph shows the diversity within each plot. Plot 2 high the highest diversity at 1.2366849, and plot 60 had the lowest diversity of 0.0836497.

spatial_div<-ggplot(data= diversity, mapping=aes(x=Plot_Number, y=div, fill=div))+
  geom_bar(stat = "identity")+
  scale_fill_viridis_c(option="turbo")+
  labs(x = "Plot Number", y = "Diversity", title ="Species Diversity Across Plots")
spatial_div

Create a Merged Dataset

Merged <- merge(canopy, plots) %>%
  merge(DBHsAndDensity)%>%
  merge(diversity) %>% 
  st_as_sf(coords = c(x="lon", y="lat")) %>%
  st_set_crs(.,value=4326) %>%
  st_transform(.,crs=st_crs(DEM))

Merged_With_TreeHeightsandAge <- merge(canopy, plots) %>%
  merge(TreeHeightsandAge) %>%
  st_as_sf(coords = c(x="lon", y="lat")) %>%
  st_set_crs(.,value=4326) %>%
  st_transform(.,crs=st_crs(DEM))

Plotted Species Diversity

Plotted_div<-ggplot(data=Merged, mapping=aes(x=plots$lon, y=plots$lat, color =div))+
  geom_point()+
  scale_color_viridis_c(option="turbo")+
   labs(x = "Longitude", y = "Latitude", title ="Species Occurene Throughout Plots", color ="Diversity")

Plotted_div

Stack All Topographic Variables

all_topo <- raster::stack(DEM,
              TRI,
              slope, aspect, 
              eastwest, northsouth,
              tpi10m, tpi100m, tpi1000m)
            
plot(all_topo)

Combine and Clean all Data at 10 Meter Resolution

envi <- raster::extract(all_topo, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean, file="proj_data/envi.tif", overwrite=TRUE
                        ) %>%     
  scale() #this makes the regression coefficients comparable 


envi<-as.data.frame(envi)

bind<- bind_cols(Merged,envi) 


#LULC Introduction 
envi2 <- raster::extract(lulc, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"), small=TRUE, fun=mean, file="proj_data/envi2.tif", overwrite=TRUE
                         ) %>%
  as.data.frame() %>%
  left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
  dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)

bighornData <- bind_cols(bind,as.data.frame(envi2)) %>%
  mutate(lulc=as.factor(lulc)) %>%
  mutate(lulc_reduced= as.factor(case_when(
    lulc==149~149,
    lulc==151~151,
    lulc==155~155,
   TRUE~1,
  )))
#  dplyr::select(-"X")

Statistics Per Plot

This image displays the data collected at each plot in a spatial format. The graphic in the top left shows the plots, numbered. The next image shows the mean Gap Fraction per plot. Followed by the mean PAR LAI per plot, the mean PAR Average per plot and the mean Canopy Density per plot. The second row of images shows the mean Leaf Angle per plot, the standard deviation of the Gap Fraction per plot, the standard deviation of the PAR LAI per plot, the standard deviation of the PAR Average per plot and the standard deviation of the Canopy Density per plot.

plot(bind) 
 title(main = "Extracting Points", line= 8)

Not Scaled- For Individual Maps!

enviNS <- raster::extract(all_topo, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean,  file="proj_data/enviNS.tif", overwrite=TRUE)

enviNS<-as.data.frame(enviNS)

bindNS<- bind_cols(Merged,enviNS) 

treesandbindNS<- merge(trees, bindNS)

#Add LULC
envi2NS <- raster::extract(lulc, st_buffer(Merged, dist = 5, endCapStyle="SQUARE"), small=TRUE, fun=mean, file="proj_data/envi2NS.tif", overwrite=TRUE
                         ) %>%
  as.data.frame() %>%
  left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
  dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)

bighornDataNS <- bind_cols(bindNS,as.data.frame(envi2NS)) %>%
  mutate(lulc=as.factor(lulc)) %>%
  mutate(lulc_reduced= as.factor(case_when(
    lulc==149~149,
    lulc==151~151,
    lulc==155~155,
    TRUE~1,
  )))

Tree Height Combine and Clean at 10 Meter Resolution

enviTree <- raster::extract(all_topo, st_buffer(Merged_With_TreeHeightsandAge,
                                                dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean, 
                            file="proj_data/enviTree.tif", overwrite=TRUE
                            )  %>%   
  scale() #this makes the regression coefficients comparable 
enviTree<-as.data.frame(enviTree) 


bindTree<- dplyr::bind_cols(Merged_With_TreeHeightsandAge,enviTree) %>%
  dplyr::rename(Average_Stand_Age=Average_Stamd_Age)

                

# LULC
envi2Tree <- raster::extract(lulc, st_buffer(Merged_With_TreeHeightsandAge, dist = 5, endCapStyle="SQUARE"), 
                             small=TRUE, fun=mean,file="proj_data/envi2Tree.tif", overwrite=TRUE) %>%
                              as.data.frame() %>%
  left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
  dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)

bighornDataTreeHeightsAndAge <- bind_cols(bindTree,as.data.frame(envi2Tree)) %>%
  mutate(lulc=as.factor(lulc)) %>%
  dplyr::select(-"X") %>%
  mutate(lulc_reduced= as.factor(case_when(
    lulc==149~149,
    lulc==151~151,
    lulc==155~155,
    TRUE~1,
  )))

Tree Height Not Scaled for Graphing

enviTreeNS <- raster::extract(all_topo, st_buffer(Merged_With_TreeHeightsandAge,
                                                dist = 5, endCapStyle="SQUARE"),small=TRUE, fun=mean, 
                            file="proj_data/enviTreeNS.tif", overwrite=TRUE
                            )  

enviTreeNS<-as.data.frame(enviTreeNS) 

bindTreeNS<- dplyr::bind_cols(Merged_With_TreeHeightsandAge,enviTreeNS) %>%
  dplyr::rename(Average_Stand_Age=Average_Stamd_Age)


#LULC
envi2TreeNS <- raster::extract(lulc, st_buffer(Merged_With_TreeHeightsandAge, dist = 5, endCapStyle="SQUARE"), small=TRUE, fun=mean,
                               file="proj_data/envi2TreeNS.tif", overwrite=TRUE) %>%
                              as.data.frame() %>%
  left_join(read_csv("proj_data/LULC.csv"), by=c("V1"="Value")) %>%
  dplyr::select(lulc=V1, LULC_Name =ECOLSYS_LU...2)

NSbighornDataTreeHeightsAndAge <- bind_cols(bindTreeNS,as.data.frame(envi2TreeNS)) %>%
  mutate(lulc=as.factor(lulc)) %>%
  dplyr::select(-"X") %>%
  mutate(lulc_reduced= as.factor(case_when(
    lulc==149~149,
    lulc==151~151,
    lulc==155~155,
    TRUE~1,
  )))

Interactive Study Site Plots

# Interactive Map of All Plots 
tmap::tmap_mode(mode='view')
tm_shape(bindTree)+
  tm_dots(col='black')
#Isolate Interactive Lidar Plots 
LidarplotsInt <-bindTree[-c(7,12,13, 23,24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 43, 51,57, 58,59, 65,66, 69),]

#Interactive Map of Lidar Plots
tmap::tmap_mode(mode='view')
tm_shape(LidarplotsInt, )+
  tm_dots(col='blue')

Forest Structure Correlations at 10 Meters

view(bighornData)

cor_bighorn10m <- bighornData %>% data.frame() %>%
  dplyr::select(
#"Plot_Number"      ,
"Gap_Fraction_mean",
"Stand_Density",
"Average_DBH",
"DBH_Sd",
"PAR_LAI_mean"     ,
"PAR_Average_mean"  ,
"Canopy_Density_mean",
"Leaf_Angle_mean" ,
"Gap_Fraction_sd",
"PAR_LAI_sd"  ,
    "PAR_Average_sd"    ,
    "Canopy_Density_sd",
    "Leaf_Angle_sd"   ,
    "div",
    "DEM"     ,
    "TRI"          ,
    "slope"        ,
    "aspect"  ,
    "eastwest"    ,
    "northsouth"       ,
    "tpi10m"          ,
    "tpi100m"          ,
    "tpi1000m"         , 
  )  %>%
  na.omit()

my_cor10<-cor(cor_bighorn10m)
view(my_cor10)
my_pmat10<-cor_pmat(cor_bighorn10m)
view(my_pmat10%>% format(scientific=F))
ggcorrplot(corr = my_cor10, p.mat = my_pmat10)

my_cor10kendall<-cor(cor_bighorn10m, method="kendall")
view(my_cor10kendall)
my_pmat10kendall<-cor_pmat(cor_bighorn10m, method="kendall")
view(my_pmat10kendall %>% format(scientific=F))
ggcorrplot(corr = my_cor10kendall, p.mat = my_pmat10kendall)

Tree Height Correlations at 10 Meeters

cor_bighornTreeHeightsAndAge <-bighornDataTreeHeightsAndAge %>% 
  data.frame() %>%
  dplyr::select(
"Mean_Tree_Height" ,
  "Sd_Tree_Height",
   "Max_Tree_Height",
    "Min_Tree_Height",
   "Average_Stand_Age",
    "DEM"     ,
    "TRI"          ,
    "slope"        ,
    "aspect"  ,
    "eastwest"    ,
    "northsouth"       ,
    "tpi10m"          ,
    "tpi100m"          ,
    "tpi1000m"         ,
  ) %>% 

  na.omit()

my_corTrees<-cor(cor_bighornTreeHeightsAndAge)
my_pmatTrees<-cor_pmat(cor_bighornTreeHeightsAndAge)
ggcorrplot(corr = my_corTrees,p.mat = my_pmatTrees)

my_corTreesKendall<-cor(cor_bighornTreeHeightsAndAge, method="kendall")
my_pmatTreesKendall<-cor_pmat(cor_bighornTreeHeightsAndAge, method="kendall")
ggcorrplot(corr = my_corTreesKendall,p.mat = my_pmatTreesKendall)

Model Building

####Gap Fraction mean ####

lm_GapFractionmMean10_exhaustive<-regsubsets(Gap_Fraction_mean~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+tpi1000m+lulc_reduced, data=bighornData,intercept=TRUE,method="exhaustive", nvmax =12,nbest =  3,)

lm_GapFractionmMean10<-lm_GapFractionmMean10_exhaustive %>% 
  broom::tidy() 
view(lm_GapFractionmMean10)


####Gap Fraction sd ####
lm_GapFraction_sd10m_exhaustive<-regsubsets(Gap_Fraction_sd~DEM+TRI+ slope+ aspect+  eastwest+ northsouth+ tpi10m+ tpi100m+tpi1000m+ lulc_reduced,data=bighornData, intercept=TRUE,  method="exhaustive", nvmax =12,  nbest =  3)

lm_GapFraction_sd10m<-lm_GapFraction_sd10m_exhaustive %>% 
  broom::tidy() 
view(lm_GapFraction_sd10m)


#####################################################################################
####PAR LAI mean ####
lm_PAR_LAI_mean_10m_exhaustive<-regsubsets(PAR_LAI_mean~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3,)

lm_PAR_LAI_mean_10m<-lm_PAR_LAI_mean_10m_exhaustive %>% 
  broom::tidy() 
view(lm_PAR_LAI_mean_10m)


####PAR LAI sd ####
lm_PAR_LAI_sd_10m_exhaustive<-regsubsets(PAR_LAI_sd~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_PAR_LAI_sd_10m<-lm_PAR_LAI_sd_10m_exhaustive %>% 
  broom::tidy() 
view(lm_PAR_LAI_sd_10m)


#####################################################################################
####PAR Average mean ####
lm_PAR_Average_mean_exhaustive<-regsubsets(PAR_Average_mean~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_PAR_Average_mean_10m<-lm_PAR_Average_mean_exhaustive %>% 
  broom::tidy() 
view(lm_PAR_Average_mean_10m)


####PAR Average SD ####
lm_PAR_Average_sd10m_exhaustive<-regsubsets(PAR_Average_sd~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_PAR_Average_sd10m<-lm_PAR_Average_sd10m_exhaustive %>% 
  broom::tidy() 
view(lm_PAR_Average_sd10m)
#####################################################################################
####Canopy Density Mean  ####
lm_Canopy_Density_mean_10m_exhaustive <-regsubsets(Canopy_Density_mean~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_Canopy_Density_mean_10m<-lm_Canopy_Density_mean_10m_exhaustive %>% 
  broom::tidy() 
view(lm_Canopy_Density_mean_10m) 

####Canopy Density SD####
lm_Canopy_Density_sd_10m_exhaustive <-regsubsets(Canopy_Density_sd~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_Canopy_Density_sd_10m<-lm_Canopy_Density_sd_10m_exhaustive %>% 
  broom::tidy() 
view(lm_Canopy_Density_sd_10m)

#####################################################################################
#### Leaf Angle mean####
lm_Leaf_Angle_mean_10_exhaustive <-regsubsets(Leaf_Angle_mean~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_Leaf_Angle_mean_10<-lm_Leaf_Angle_mean_10_exhaustive %>% 
  broom::tidy() 
view(lm_Leaf_Angle_mean_10)

####Leaf Angle sd ####
lm_Leaf_Angle_sd_10_exhaustive <-regsubsets(Leaf_Angle_sd~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_Leaf_Angle_sd_10<-lm_Leaf_Angle_sd_10_exhaustive %>% 
  broom::tidy() 
view(lm_Leaf_Angle_sd_10)

#####################################################################################
####Diversity ####

lm_Diversity_exhaustive <-regsubsets(div~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_Diversity_10m<-lm_Diversity_exhaustive %>% 
  broom::tidy() 
view(lm_Diversity_10m)
#####################################################################################
####Average DBH ####
lm_DBH_mean_10m_exhaustive <-regsubsets(Average_DBH~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_DBH_mean_10m<-lm_DBH_mean_10m_exhaustive %>% 
  broom::tidy() 
view(lm_DBH_mean_10m)

#####################################################################################
####Dbh SD ####
lm_DBH_sd_10m_exhaustive <-regsubsets(DBH_Sd~ DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_DBH_sd_10m<-lm_DBH_sd_10m_exhaustive %>% 
  broom::tidy() 
view(lm_DBH_sd_10m)
#####################################################################################
####Stand Density ####
lm_Stand_Density_10m_exhaustive <-regsubsets(Stand_Density~DEM+ TRI+ slope+ aspect+ eastwest+ northsouth+ tpi10m+ tpi100m+ tpi1000m+lulc_reduced,data=bighornData,intercept=TRUE,method="exhaustive", nvmax =10, nbest =  3)

lm_Stand_Density_10m<-lm_Stand_Density_10m_exhaustive %>% 
  broom::tidy() 
view(lm_Stand_Density_10m)

###################################################################################
####Mean Tree Height  ####
lm_TreeHeightMean_exhaustive<-regsubsets(Mean_Tree_Height~DEM+TRI+slope+aspect+ eastwest+northsouth+ tpi10m+  tpi100m+  tpi1000m+lulc_reduced, data=bighornDataTreeHeightsAndAge, intercept=TRUE, method="exhaustive", nvmax =10, nbest =  3)

lm_TreeHeightMean<-lm_TreeHeightMean_exhaustive %>% 
  broom::tidy() 
view(lm_TreeHeightMean)


###################################################################################
####SD Tree Height ####

lm_TreeHeight_SD_exhaustive<-regsubsets(Sd_Tree_Height~DEM+TRI+slope+aspect+ eastwest+northsouth+ tpi10m+  tpi100m+  tpi1000m+lulc_reduced, data=bighornDataTreeHeightsAndAge, intercept=TRUE, method="exhaustive", nvmax =10, nbest =  3)

lm_TreeHeightSD<-lm_TreeHeight_SD_exhaustive %>% 
  broom::tidy() 
view(lm_TreeHeightSD)

##################################################################################
#Max Tree Height
lm_Max_TreeHeight_exhaustive<-regsubsets(Max_Tree_Height~DEM+TRI+slope+aspect+ eastwest+northsouth+ tpi10m+  tpi100m+  tpi1000m+lulc_reduced, data=bighornDataTreeHeightsAndAge, intercept=TRUE, method="exhaustive", nvmax =10, nbest =  3)

lm_Max_TreeHeight<-lm_Max_TreeHeight_exhaustive %>% 
  broom::tidy() 
view(lm_Max_TreeHeight)

##################################################################################
#Average Stand Age (tpi100 and 1000 individually, and Landforms100 and 1000)
lm_Average_Stand_Age10m_exhaustive<-regsubsets(Average_Stand_Age~DEM+TRI+slope+aspect+ eastwest+northsouth+ tpi10m+  tpi100m+  tpi1000m+lulc_reduced, data=bighornDataTreeHeightsAndAge, intercept=TRUE, method="exhaustive", nvmax =10, nbest =  3)

lm_Average_Stand_Age10m<-lm_Average_Stand_Age10m_exhaustive %>% 
  broom::tidy() 
view(lm_Average_Stand_Age10m)

Forest Structure Regressions at 10 Meters

#Gap Fraction mean 
#p-value: 2.378e-05
GapFraction_mean_10m<-lm(Gap_Fraction_mean~ TRI+ tpi1000m+lulc_reduced,
                         data=bighornData)
summary(GapFraction_mean_10m)
## 
## Call:
## lm(formula = Gap_Fraction_mean ~ TRI + tpi1000m + lulc_reduced, 
##     data = bighornData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90669 -0.16654  0.00984  0.20381  0.72175 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.11244    0.09494  11.717  < 2e-16 ***
## TRI              0.07958    0.03974   2.003   0.0493 *  
## tpi1000m         0.08333    0.04003   2.082   0.0412 *  
## lulc_reduced149  0.49340    0.10851   4.547 2.39e-05 ***
## lulc_reduced151  0.60939    0.14439   4.220 7.61e-05 ***
## lulc_reduced155  0.16545    0.15041   1.100   0.2753    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3257 on 66 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.349,  Adjusted R-squared:  0.2997 
## F-statistic: 7.076 on 5 and 66 DF,  p-value: 2.378e-05
broom::tidy(GapFraction_mean_10m)
## # A tibble: 6 × 5
##   term            estimate std.error statistic  p.value
##   <chr>              <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)       1.11      0.0949     11.7  8.91e-18
## 2 TRI               0.0796    0.0397      2.00 4.93e- 2
## 3 tpi1000m          0.0833    0.0400      2.08 4.12e- 2
## 4 lulc_reduced149   0.493     0.109       4.55 2.39e- 5
## 5 lulc_reduced151   0.609     0.144       4.22 7.61e- 5
## 6 lulc_reduced155   0.165     0.150       1.10 2.75e- 1
broom::glance(GapFraction_mean_10m)
## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic   p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>     <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.349         0.300 0.326      7.08 0.0000238     5  -18.3  50.5  66.5
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
broom::augment(GapFraction_mean_10m)
## # A tibble: 72 × 11
##    .rownames Gap_Fraction_mean     TRI tpi1000m lulc_reduced .fitted  .resid
##    <chr>                 <dbl>   <dbl>    <dbl> <fct>          <dbl>   <dbl>
##  1 1                      1.22 -0.973    0.920  1               1.11  0.108 
##  2 2                      1.72 -1.41     0.545  149             1.54  0.181 
##  3 3                      1.72 -0.588    0.771  149             1.62  0.0967
##  4 4                      1.7  -1.08     1.78   149             1.67  0.0315
##  5 5                      1.78  2.53     1.24   149             1.91 -0.130 
##  6 6                      1.78 -0.944    0.349  149             1.56  0.220 
##  7 8                      1.72 -1.07    -0.0324 149             1.52  0.202 
##  8 9                      1.56 -0.182   -0.0256 149             1.59 -0.0292
##  9 10                     1.46  0.433    0.687  149             1.70 -0.238 
## 10 11                     1.58  0.0482   0.117  149             1.62 -0.0394
## # … with 62 more rows, and 4 more variables: .hat <dbl>, .sigma <dbl>,
## #   .cooksd <dbl>, .std.resid <dbl>
densityPlot(GapFraction_mean_10m$residuals)

AIC(GapFraction_mean_10m)
## [1] 50.52386
BIC(GapFraction_mean_10m)
## [1] 66.46053
# 
# MVGapFraction_mean_10m<-glm(Gap_Fraction_mean~TRI+ tpi1000m+lulc_reduced,
#                 family=gaussian(link="log"),
#                 data=bighornData)
# summary(MVGapFraction_mean_10m)
#BIC(GapFraction_mean_10m, MVGapFraction_mean_10m)

#Gap Fraction sd 
#p-value: 0.05405
# GapFraction_sd_10m<-glm(Gap_Fraction_sd~  DEM,
#                         family = gaussian(),
#                         data=bighornData)
# summary(GapFraction_sd_10m)
#densityPlot(GapFraction_sd_10m$residuals)

MVGapFraction_sd_10m<-glm(Gap_Fraction_sd~DEM,
                family=gaussian(link="log"),
                data=bighornData) 
summary(MVGapFraction_sd_10m)
## 
## Call:
## glm(formula = Gap_Fraction_sd ~ DEM, family = gaussian(link = "log"), 
##     data = bighornData)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.22114  -0.13160  -0.06747   0.07563   0.78255  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.24112    0.07891 -15.728   <2e-16 ***
## DEM         -0.14521    0.06412  -2.265   0.0266 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.03671833)
## 
##     Null deviance: 2.7575  on 72  degrees of freedom
## Residual deviance: 2.6070  on 71  degrees of freedom
## AIC: -30.09
## 
## Number of Fisher Scoring iterations: 6
AIC(MVGapFraction_sd_10m)
## [1] -30.09
BIC(MVGapFraction_sd_10m)
## [1] -23.21862
broom::tidy(MVGapFraction_sd_10m) 
## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   -1.24     0.0789    -15.7  7.82e-25
## 2 DEM           -0.145    0.0641     -2.26 2.66e- 2
broom::glance(MVGapFraction_sd_10m)
## # A tibble: 1 × 8
##   null.deviance df.null logLik   AIC   BIC deviance df.residual  nobs
##           <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int> <int>
## 1          2.76      72   18.0 -30.1 -23.2     2.61          71    73
densityPlot(MVGapFraction_sd_10m$residuals)

#####################################################################################
####PAR LAI mean ####
# #p-value: 0.1319
# PAR_LAI_mean_10m<-lm(PAR_LAI_mean~ DEM+tpi10m+lulc_reduced,
#                        data=bighornData)
# summary(PAR_LAI_mean_10m)

mvPAR_LAI_mean_10m<-glm(PAR_LAI_mean~ DEM+tpi10m+lulc_reduced,
                family=gaussian(link="log"),
                data=bighornData)
summary(mvPAR_LAI_mean_10m)
## 
## Call:
## glm(formula = PAR_LAI_mean ~ DEM + tpi10m + lulc_reduced, family = gaussian(link = "log"), 
##     data = bighornData)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4816  -0.9334  -0.2222   0.9833   4.3572  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.72732    0.20962   3.470 0.000916 ***
## DEM              0.17283    0.08178   2.113 0.038297 *  
## tpi10m           0.12018    0.08230   1.460 0.148864    
## lulc_reduced149  0.22686    0.22516   1.008 0.317289    
## lulc_reduced151  0.16763    0.27760   0.604 0.547980    
## lulc_reduced155  0.60219    0.26272   2.292 0.025045 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 2.234312)
## 
##     Null deviance: 172.23  on 72  degrees of freedom
## Residual deviance: 149.70  on 67  degrees of freedom
## AIC: 273.59
## 
## Number of Fisher Scoring iterations: 8
densityPlot(mvPAR_LAI_mean_10m$residuals)

AIC(mvPAR_LAI_mean_10m)
## [1] 273.5912
BIC(mvPAR_LAI_mean_10m) 
## [1] 289.6244
#                     df      BIC
# PAR_LAI_mean_10m    7 290.8233
# mvPAR_LAI_mean_10m  7 289.6244

####PAR LAI sd ####
#p-value: 0.02881 
PAR_LAI_sd_10m<-lm(PAR_LAI_sd~ eastwest,
                   data=bighornData)
summary(PAR_LAI_sd_10m)
## 
## Call:
## lm(formula = PAR_LAI_sd ~ eastwest, data = bighornData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73833 -0.27252 -0.06329  0.20379  1.57294 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.64220    0.04844  13.259   <2e-16 ***
## eastwest    -0.10883    0.04877  -2.231   0.0288 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4138 on 71 degrees of freedom
## Multiple R-squared:  0.06553,    Adjusted R-squared:  0.05237 
## F-statistic: 4.979 on 1 and 71 DF,  p-value: 0.02881
AIC(PAR_LAI_sd_10m)
## [1] 82.32141
BIC(PAR_LAI_sd_10m) 
## [1] 89.19279
# mvPAR_LAI_sd_10m<-glm(PAR_LAI_sd~ eastwest,
#                 family=gaussian(link="log"),
#                 data=bighornData)
# summary(mvPAR_LAI_sd10m)
 # BIC(PAR_LAI_sd_10m,mvPAR_LAI_sd_10m) 
#                  df      BIC
# PAR_LAI_sd_10m    3 89.19279
# mvPAR_LAI_sd_10m  3 89.71298


#####################################################################################
####PAR Average mean ####
#p-value: 0.008254
PAR_Average_mean_10m<-lm(PAR_Average_mean~  lulc_reduced,
                         data=bighornData)
summary(PAR_Average_mean_10m)
## 
## Call:
## lm(formula = PAR_Average_mean ~ lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -410.28 -127.40  -67.51  107.00  777.92 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       420.48      61.84   6.799 3.06e-09 ***
## lulc_reduced149  -227.87      69.76  -3.266   0.0017 ** 
## lulc_reduced151  -286.57      94.46  -3.034   0.0034 ** 
## lulc_reduced155  -219.35      97.78  -2.243   0.0281 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 214.2 on 69 degrees of freedom
## Multiple R-squared:  0.1557, Adjusted R-squared:  0.119 
## F-statistic: 4.241 on 3 and 69 DF,  p-value: 0.008254
densityPlot(PAR_Average_mean_10m$residuals)

AIC(PAR_Average_mean_10m)
## [1] 996.6328
BIC(PAR_Average_mean_10m) 
## [1] 1008.085
# mvPAR_Average_mean_10m<-glm(PAR_Average_mean~  lulc_reduced,
#                 family=gaussian(link="log"),
#                 data=bighornData)
# summary(mvPAR_Average_mean_10m)
# BIC(PAR_Average_mean_10m,mvPAR_Average_mean_10m)               
#                         df      BIC
# PAR_Average_mean_10m    5 1008.085
# mvPAR_Average_mean_10m  5 1008.085


####PAR Average SD ####
#p-value: 0.01623
# PAR_Average_sd10m<-lm(PAR_Average_sd~ TRI+slope+ northsouth,
#                    data=bighornData)
# summary(PAR_Average_sd10m)

mvPAR_Average_sd10m<-glm(PAR_Average_sd~ TRI+slope+ northsouth,
                family=gaussian(link="log"),
                data=bighornData)
summary(mvPAR_Average_sd10m)
## 
## Call:
## glm(formula = PAR_Average_sd ~ TRI + slope + northsouth, family = gaussian(link = "log"), 
##     data = bighornData)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -262.81  -104.18    -8.98    96.41   539.49  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.7304     0.2414  19.593   <2e-16 ***
## TRI           1.6649     0.8282   2.010   0.0483 *  
## slope        -2.2299     0.9365  -2.381   0.0200 *  
## northsouth   -0.3935     0.1950  -2.018   0.0475 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 28999.74)
## 
##     Null deviance: 2451433  on 72  degrees of freedom
## Residual deviance: 2000577  on 69  degrees of freedom
## AIC: 963.11
## 
## Number of Fisher Scoring iterations: 17
AIC(mvPAR_Average_sd10m)
## [1] 963.1146
BIC(mvPAR_Average_sd10m) 
## [1] 974.5668
# BIC(PAR_Average_sd10m,mvPAR_Average_sd10m)
#                     df      BIC
# PAR_Average_sd10m    5 978.5905
# mvPAR_Average_sd10m  5 974.5668

#####################################################################################
#### Canopy Density Mean ####
#p-value: 8.287e-06
Canopy_Density_mean_10m<-lm(Canopy_Density_mean~  DEM+ lulc_reduced,
                      data=bighornData)
summary(Canopy_Density_mean_10m)
## 
## Call:
## lm(formula = Canopy_Density_mean ~ DEM + lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.178  -2.206   1.187   4.191  22.762 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       55.305      2.925  18.908  < 2e-16 ***
## DEM                3.612      1.328   2.721  0.00826 ** 
## lulc_reduced149   14.206      3.317   4.283 5.93e-05 ***
## lulc_reduced151   12.912      4.605   2.804  0.00658 ** 
## lulc_reduced155    8.106      4.715   1.719  0.09010 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.06 on 68 degrees of freedom
## Multiple R-squared:  0.3421, Adjusted R-squared:  0.3034 
## F-statistic: 8.841 on 4 and 68 DF,  p-value: 8.287e-06
AIC(Canopy_Density_mean_10m)
## [1] 551.0652
BIC(Canopy_Density_mean_10m)
## [1] 564.8079
densityPlot(Canopy_Density_mean_10m$residuals)

# mvCanopy_Density_mean_10m<-glm(Canopy_Density_mean~  DEM+ lulc_reduced,
#                 family=gaussian(link="logit"),
#                 data=bighornData)
# summary(mvCanopy_Density_mean_10m)
#BIC(Canopy_Density_mean_10m,mvCanopy_Density_mean_10m) 
#                           df      BIC
# Canopy_Density_mean_10m    6 564.8079
# mvCanopy_Density_mean_10m  6 566.1761

#####Canopy Density SD ####
#5.175e-05
# Canopy_Density_sd_10m<-lm(Canopy_Density_sd~ DEM+tpi100m+ lulc_reduced,
#                             data=bighornData)
# summary(Canopy_Density_sd_10m)

mvCanopy_Density_sd_10m<-glm(Canopy_Density_sd~ DEM+tpi100m+ lulc_reduced,
                family=gaussian(link="log"),
                data=bighornData)
summary(mvCanopy_Density_sd_10m)
## 
## Call:
## glm(formula = Canopy_Density_sd ~ DEM + tpi100m + lulc_reduced, 
##     family = gaussian(link = "log"), data = bighornData)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -7.4643  -2.3473  -1.1044   0.4404  14.6721  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.23457    0.15057  14.841  < 2e-16 ***
## DEM             -0.16665    0.07603  -2.192  0.03186 *  
## tpi100m         -0.25758    0.09212  -2.796  0.00674 ** 
## lulc_reduced149 -0.74751    0.22374  -3.341  0.00137 ** 
## lulc_reduced151 -0.88263    0.45086  -1.958  0.05444 .  
## lulc_reduced155 -0.61829    0.28028  -2.206  0.03082 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 22.61968)
## 
##     Null deviance: 2338.1  on 72  degrees of freedom
## Residual deviance: 1515.5  on 67  degrees of freedom
## AIC: 442.58
## 
## Number of Fisher Scoring iterations: 8
AIC(mvCanopy_Density_sd_10m)
## [1] 442.5765
BIC(mvCanopy_Density_sd_10m)
## [1] 458.6098
#BIC(Canopy_Density_sd_10m,mvCanopy_Density_sd_10m) 
#                         df      BIC
# Canopy_Density_sd_10m    7 461.2300
# mvCanopy_Density_sd_10m  7 458.6098

####################################################################################
#### Leaf Angle mean ####
#p-value: 0.5935
Leaf_Angle_mean_10<-lm(Leaf_Angle_mean~ lulc_reduced,
                          data=bighornData)
summary(Leaf_Angle_mean_10)
## 
## Call:
## lm(formula = Leaf_Angle_mean ~ lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.627  -8.346   1.694   9.114  35.162 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       48.727      4.390  11.100   <2e-16 ***
## lulc_reduced149    5.660      4.952   1.143    0.257    
## lulc_reduced151    1.049      6.706   0.156    0.876    
## lulc_reduced155    1.311      6.941   0.189    0.851    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.21 on 69 degrees of freedom
## Multiple R-squared:  0.02696,    Adjusted R-squared:  -0.01534 
## F-statistic: 0.6374 on 3 and 69 DF,  p-value: 0.5935
AIC(Leaf_Angle_mean_10)
## [1] 610.4254
BIC(Leaf_Angle_mean_10)
## [1] 621.8777
densityPlot(Leaf_Angle_mean_10$residuals)

#### Leaf Angle sd ####
#p-value: 0.001394
Leaf_Angle_sd_10<-lm(Leaf_Angle_sd~ tpi1000m+lulc_reduced,
                       data=bighornData)
summary(Leaf_Angle_sd_10)
## 
## Call:
## lm(formula = Leaf_Angle_sd ~ tpi1000m + lulc_reduced, data = bighornData)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.8374  -6.5164  -0.5118   7.8988  24.1936 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      24.7175     2.9241   8.453 3.67e-12 ***
## tpi1000m         -4.1758     1.2449  -3.354  0.00131 ** 
## lulc_reduced149  -6.5279     3.3152  -1.969  0.05308 .  
## lulc_reduced151  -0.9586     4.4813  -0.214  0.83126    
## lulc_reduced155  -3.1830     4.6759  -0.681  0.49839    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.13 on 67 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2296, Adjusted R-squared:  0.1836 
## F-statistic: 4.991 on 4 and 67 DF,  p-value: 0.001394
AIC(Leaf_Angle_sd_10)
## [1] 544.5641
BIC(Leaf_Angle_sd_10)
## [1] 558.2241
densityPlot(Leaf_Angle_sd_10$residuals)

#####################################################################################
####Diversity ####
# p-value: 0.0009821
Diversity_10m<-lm(div~DEM+northsouth+lulc_reduced,
                     data=bighornData)
summary(Diversity_10m) 
## 
## Call:
## lm(formula = div ~ DEM + northsouth + lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3811 -0.1322 -0.0067  0.1049  0.7027 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.38439    0.06148   6.252 3.22e-08 ***
## DEM              0.09229    0.02705   3.412   0.0011 ** 
## northsouth      -0.04444    0.02512  -1.769   0.0815 .  
## lulc_reduced149  0.16987    0.06971   2.437   0.0175 *  
## lulc_reduced151  0.03923    0.09604   0.408   0.6842    
## lulc_reduced155  0.15174    0.09824   1.545   0.1272    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.205 on 67 degrees of freedom
## Multiple R-squared:  0.2594, Adjusted R-squared:  0.2041 
## F-statistic: 4.694 on 5 and 67 DF,  p-value: 0.0009821
AIC(Diversity_10m)
## [1] -16.45697
BIC(Diversity_10m)
## [1] -0.4237531
densityPlot(Diversity_10m$residuals)

#####################################################################################

#### Average DBH ####
#p-value: 0.0009843
# DBH_mean_10m<-lm(Average_DBH~ tpi100m+ lulc_reduced,
#                   data=bighornData)
# summary(DBH_mean_10m)

mvDBH_mean_10m<-glm(Average_DBH~ tpi100m+ lulc_reduced,
                family=gaussian(link="log"),
                data=bighornData)
summary(mvDBH_mean_10m)
## 
## Call:
## glm(formula = Average_DBH ~ tpi100m + lulc_reduced, family = gaussian(link = "log"), 
##     data = bighornData)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -18.515   -5.417   -1.854    5.174   28.218  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.44753    0.08187  42.110  < 2e-16 ***
## tpi100m         -0.11882    0.03845  -3.091  0.00289 ** 
## lulc_reduced149 -0.25190    0.09853  -2.557  0.01281 *  
## lulc_reduced151  0.10019    0.11714   0.855  0.39542    
## lulc_reduced155 -0.29661    0.15219  -1.949  0.05542 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 77.95788)
## 
##     Null deviance: 7049.7  on 72  degrees of freedom
## Residual deviance: 5301.1  on 68  degrees of freedom
## AIC: 531.99
## 
## Number of Fisher Scoring iterations: 5
AIC(mvDBH_mean_10m)
## [1] 531.9853
BIC(mvDBH_mean_10m)
## [1] 545.7281
#BIC(DBH_mean_10m,mvDBH_mean_10m)
#                df      BIC
# DBH_mean_10m    6 546.9549
# mvDBH_mean_10m  6 545.7281

#####################################################################################

#### SD DBH #### 
#p-value: 0.001216
# DBH_sd_10m<-lm(DBH_Sd~ tpi1000m+ lulc_reduced,
#                  data=bighornData)
# summary(DBH_sd_10m)


mvDBH_sd_10m<-glm(DBH_Sd~ tpi1000m+ lulc_reduced,
                family=gaussian(link="log"), 
                start = c(1,1,1,1,1),
                data=bighornData)
summary(mvDBH_sd_10m)
## 
## Call:
## glm(formula = DBH_Sd ~ tpi1000m + lulc_reduced, family = gaussian(link = "log"), 
##     data = bighornData, start = c(1, 1, 1, 1, 1))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -10.9495   -2.9485   -0.4068    3.0472   15.9326  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1.91597    0.22261   8.607 1.94e-12 ***
## tpi1000m        -0.20102    0.07326  -2.744  0.00778 ** 
## lulc_reduced149  0.14843    0.24541   0.605  0.54732    
## lulc_reduced151  0.76101    0.24972   3.047  0.00330 ** 
## lulc_reduced155 -0.22064    0.37252  -0.592  0.55565    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 29.08617)
## 
##     Null deviance: 2658.4  on 71  degrees of freedom
## Residual deviance: 1948.7  on 67  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 453.8
## 
## Number of Fisher Scoring iterations: 10
AIC(mvDBH_sd_10m)
## [1] 453.803
BIC(mvDBH_sd_10m)
## [1] 467.463
densityPlot(mvDBH_sd_10m$residuals)

#BIC(DBH_sd_10m,mvDBH_sd_10m)
#              df      BIC
# DBH_sd_10m    6 470.7204
# mvDBH_sd_10m  6 467.4630

#####################################################################################
#S####tand Density ####
#p-value: 1.806e-06
Stand_Density_10m<-lm(Stand_Density~ DEM+ tpi100m+ lulc_reduced,
               data=bighornData)
summary(Stand_Density_10m) 
## 
## Call:
## lm(formula = Stand_Density ~ DEM + tpi100m + lulc_reduced, data = bighornData)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9979 -1.8267 -0.1134  1.6456  9.6533 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.0015     0.9799   6.125 5.38e-08 ***
## DEM               1.2977     0.4563   2.844  0.00591 ** 
## tpi100m           1.3722     0.4234   3.241  0.00186 ** 
## lulc_reduced149   3.0350     1.1082   2.739  0.00790 ** 
## lulc_reduced151  -0.2019     1.5491  -0.130  0.89671    
## lulc_reduced155   1.5211     1.5988   0.951  0.34483    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.36 on 67 degrees of freedom
## Multiple R-squared:  0.3969, Adjusted R-squared:  0.3519 
## F-statistic: 8.819 on 5 and 67 DF,  p-value: 1.806e-06
AIC(Stand_Density_10m)
## [1] 391.8308
BIC(Stand_Density_10m)
## [1] 407.864
densityPlot(Stand_Density_10m$residuals)

#####################################################################################

Tree Height Regressions at 10 Meters

#Mean Tree Height
#p-value: 0.0199
 Mean_Tree_Height10m<-lm(Mean_Tree_Height~northsouth+lulc_reduced,
                     data=bighornDataTreeHeightsAndAge)
summary(Mean_Tree_Height10m)
## 
## Call:
## lm(formula = Mean_Tree_Height ~ northsouth + lulc_reduced, data = bighornDataTreeHeightsAndAge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.9490 -2.0172 -0.4666  2.3026  9.2168 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       7.3446     1.3570   5.412 3.64e-06 ***
## northsouth       -0.8025     0.5636  -1.424  0.16268    
## lulc_reduced149   5.3506     1.5301   3.497  0.00122 ** 
## lulc_reduced151   5.5459     2.3165   2.394  0.02170 *  
## lulc_reduced155   6.9510     3.6681   1.895  0.06572 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.43 on 38 degrees of freedom
##   (30 observations deleted due to missingness)
## Multiple R-squared:  0.259,  Adjusted R-squared:  0.181 
## F-statistic: 3.321 on 4 and 38 DF,  p-value: 0.0199
AIC(Mean_Tree_Height10m)
## [1] 234.7079
BIC(Mean_Tree_Height10m)
## [1] 245.2751
densityPlot( Mean_Tree_Height10m$residuals)

###################################################################################
#SD Tree Height
#p-value: 0.006275
SD_Tree_Height10m<-lm(Sd_Tree_Height~ tpi100m+lulc_reduced,
                      data=bighornDataTreeHeightsAndAge)
summary(SD_Tree_Height10m)
## 
## Call:
## lm(formula = Sd_Tree_Height ~ tpi100m + lulc_reduced, data = bighornDataTreeHeightsAndAge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0817 -1.2138 -0.1103  0.5723  4.4815 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       5.0266     0.6550   7.674 3.07e-09 ***
## tpi100m          -0.8031     0.2634  -3.049  0.00417 ** 
## lulc_reduced149  -1.0080     0.7271  -1.386  0.17371    
## lulc_reduced151   0.7790     1.0901   0.715  0.47922    
## lulc_reduced155  -0.8551     1.8526  -0.462  0.64704    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.732 on 38 degrees of freedom
##   (30 observations deleted due to missingness)
## Multiple R-squared:  0.308,  Adjusted R-squared:  0.2352 
## F-statistic: 4.229 on 4 and 38 DF,  p-value: 0.006275
AIC(SD_Tree_Height10m)
## [1] 175.9737
BIC(SD_Tree_Height10m)
## [1] 186.5409
densityPlot(SD_Tree_Height10m$residuals)

##################################################################################
#Max Tree Height (northsouth aspect)
#p-value: 0.00969
Max_Tree_Height10m<-lm(Max_Tree_Height~ aspect+ tpi1000m +lulc_reduced,
                       data=bighornDataTreeHeightsAndAge)
summary(Max_Tree_Height10m)
## 
## Call:
## lm(formula = Max_Tree_Height ~ aspect + tpi1000m + lulc_reduced, 
##     data = bighornDataTreeHeightsAndAge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6333 -2.8068  0.3041  1.9522  6.8380 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      15.8762     1.3137  12.085 2.07e-14 ***
## aspect            0.8112     0.5026   1.614   0.1150    
## tpi1000m         -1.1283     0.5897  -1.913   0.0635 .  
## lulc_reduced149   2.3294     1.4470   1.610   0.1159    
## lulc_reduced151   5.5133     2.1752   2.535   0.0156 *  
## lulc_reduced155   6.5971     3.6305   1.817   0.0773 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.388 on 37 degrees of freedom
##   (30 observations deleted due to missingness)
## Multiple R-squared:  0.326,  Adjusted R-squared:  0.235 
## F-statistic:  3.58 on 5 and 37 DF,  p-value: 0.00969
AIC(Max_Tree_Height10m)
## [1] 234.5193
BIC(Max_Tree_Height10m)
## [1] 246.8477
densityPlot(Max_Tree_Height10m$residuals)

##################################################################################
#Average Stand Age 
#p-value: 0.0001305
Average_Stand_Age10m<-lm(Average_Stand_Age~aspect+ tpi1000m +lulc_reduced,
                         data=bighornDataTreeHeightsAndAge)
summary(Average_Stand_Age10m)
## 
## Call:
## lm(formula = Average_Stand_Age ~ aspect + tpi1000m + lulc_reduced, 
##     data = bighornDataTreeHeightsAndAge)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.865 -10.901  -2.243   8.560  56.634 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       74.491      5.781  12.885  < 2e-16 ***
## aspect             5.211      2.422   2.151  0.03510 *  
## tpi1000m          -6.481      2.449  -2.646  0.01016 *  
## lulc_reduced149  -19.671      6.534  -3.010  0.00369 ** 
## lulc_reduced151    1.676      8.864   0.189  0.85066    
## lulc_reduced155  -18.959      9.432  -2.010  0.04852 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.89 on 66 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.3113, Adjusted R-squared:  0.2591 
## F-statistic: 5.965 on 5 and 66 DF,  p-value: 0.0001305
AIC(Average_Stand_Age10m)
## [1] 642.6376
BIC(Average_Stand_Age10m)
## [1] 658.5743
densityPlot(Average_Stand_Age10m$residuals)

Regression Results

Regression Results

#Gap Fraction Mean 
par(mfrow=c(2,3))
plot(GapFraction_mean_10m, 1:6)

gvlma(GapFraction_mean_10m)
## 
## Call:
## lm(formula = Gap_Fraction_mean ~ TRI + tpi1000m + lulc_reduced, 
##     data = bighornData)
## 
## Coefficients:
##     (Intercept)              TRI         tpi1000m  lulc_reduced149  
##         1.11244          0.07958          0.08333          0.49340  
## lulc_reduced151  lulc_reduced155  
##         0.60939          0.16545  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = GapFraction_mean_10m) 
## 
##                     Value p-value                   Decision
## Global Stat        12.885 0.01185 Assumptions NOT satisfied!
## Skewness            1.457 0.22742    Assumptions acceptable.
## Kurtosis            1.715 0.19028    Assumptions acceptable.
## Link Function       4.881 0.02716 Assumptions NOT satisfied!
## Heteroscedasticity  4.832 0.02794 Assumptions NOT satisfied!
# Coefficients:
#     (Intercept)              TRI         tpi1000m  lulc_reduced149 lulc_reduced151  lulc_reduced155    
#         1.11244          0.07958          0.08333          0.49340          0.60939          0.16545  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = GapFraction_mean_10m) 
# 
#                     Value p-value                   Decision
# Global Stat        12.885 0.01185 Assumptions NOT satisfied!
# Skewness            1.457 0.22742    Assumptions acceptable.
# Kurtosis            1.715 0.19028    Assumptions acceptable.
# Link Function       4.881 0.02716 Assumptions NOT satisfied!
# Heteroscedasticity  4.832 0.02794 Assumptions NOT satisfied!

#Gap Fraction SD
par(mfrow=c(2,3))
plot(MVGapFraction_sd_10m, 1:6)

#PAR LAI Mean
par(mfrow=c(2,3))
plot(mvPAR_LAI_mean_10m, 1:6)

#PAR LAI SD
par(mfrow=c(2,3))
plot(PAR_LAI_sd_10m, 1:6) 

gvlma(PAR_LAI_sd_10m)
## 
## Call:
## lm(formula = PAR_LAI_sd ~ eastwest, data = bighornData)
## 
## Coefficients:
## (Intercept)     eastwest  
##      0.6422      -0.1088  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = PAR_LAI_sd_10m) 
## 
##                      Value   p-value                   Decision
## Global Stat        31.2255 2.754e-06 Assumptions NOT satisfied!
## Skewness           15.6306 7.700e-05 Assumptions NOT satisfied!
## Kurtosis           14.2828 1.573e-04 Assumptions NOT satisfied!
## Link Function       0.7860 3.753e-01    Assumptions acceptable.
## Heteroscedasticity  0.5261 4.683e-01    Assumptions acceptable.
# Call:
# lm(formula = PAR_LAI_sd ~ eastwest, data = bighornData)
# 
# Coefficients:
# (Intercept)     eastwest  
#      0.6422      -0.1088  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = PAR_LAI_sd_10m) 
# 
#                      Value   p-value                   Decision
# Global Stat        31.2255 2.754e-06 Assumptions NOT satisfied!
# Skewness           15.6306 7.700e-05 Assumptions NOT satisfied!
# Kurtosis           14.2828 1.573e-04 Assumptions NOT satisfied!
# Link Function       0.7860 3.753e-01    Assumptions acceptable.
# Heteroscedasticity  0.5261 4.683e-01    Assumptions acceptable.


#PAR Average Mean
par(mfrow=c(2,3))
plot(PAR_Average_mean_10m, 1:6)

gvlma(PAR_Average_mean_10m)
## 
## Call:
## lm(formula = PAR_Average_mean ~ lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)  lulc_reduced149  lulc_reduced151  lulc_reduced155  
##           420.5           -227.9           -286.6           -219.4  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = PAR_Average_mean_10m) 
## 
##                        Value   p-value                   Decision
## Global Stat        4.265e+01 1.222e-08 Assumptions NOT satisfied!
## Skewness           1.811e+01 2.089e-05 Assumptions NOT satisfied!
## Kurtosis           1.321e+01 2.779e-04 Assumptions NOT satisfied!
## Link Function      1.853e-14 1.000e+00    Assumptions acceptable.
## Heteroscedasticity 1.133e+01 7.623e-04 Assumptions NOT satisfied!
# Call:
# lm(formula = PAR_Average_mean ~ lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)  lulc_reduced149  lulc_reduced151  lulc_reduced155  
#           420.5           -227.9           -286.6           -219.4  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = PAR_Average_mean_10m) 
# 
#                        Value   p-value                   Decision
# Global Stat        4.265e+01 1.222e-08 Assumptions NOT satisfied!
# Skewness           1.811e+01 2.089e-05 Assumptions NOT satisfied!
# Kurtosis           1.321e+01 2.779e-04 Assumptions NOT satisfied!
# Link Function      1.853e-14 1.000e+00    Assumptions acceptable.
# Heteroscedasticity 1.133e+01 7.623e-04 Assumptions NOT satisfied!

#PAR Average SD
par(mfrow=c(2,3))
plot(mvPAR_Average_sd10m, 1:6)

#Canopy Density Mean
par(mfrow=c(2,3))
plot(Canopy_Density_mean_10m, 1:6)

gvlma(Canopy_Density_mean_10m)
## 
## Call:
## lm(formula = Canopy_Density_mean ~ DEM + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)              DEM  lulc_reduced149  lulc_reduced151  
##          55.305            3.612           14.206           12.912  
## lulc_reduced155  
##           8.106  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Canopy_Density_mean_10m) 
## 
##                     Value   p-value                   Decision
## Global Stat        87.465 0.000e+00 Assumptions NOT satisfied!
## Skewness           23.079 1.555e-06 Assumptions NOT satisfied!
## Kurtosis           51.359 7.694e-13 Assumptions NOT satisfied!
## Link Function       4.812 2.826e-02 Assumptions NOT satisfied!
## Heteroscedasticity  8.216 4.153e-03 Assumptions NOT satisfied!
# Call:
# lm(formula = Canopy_Density_mean ~ DEM + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)              DEM  lulc_reduced149  lulc_reduced151  lulc_reduced155  
#          55.305            3.612           14.206           12.912            8.106  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Canopy_Density_mean_10m) 
# 
#                     Value   p-value                   Decision
# Global Stat        87.465 0.000e+00 Assumptions NOT satisfied!
# Skewness           23.079 1.555e-06 Assumptions NOT satisfied!
# Kurtosis           51.359 7.694e-13 Assumptions NOT satisfied!
# Link Function       4.812 2.826e-02 Assumptions NOT satisfied!
# Heteroscedasticity  8.216 4.153e-03 Assumptions NOT satisfied!


#Canopy Density SD
par(mfrow=c(2,3))
plot(mvCanopy_Density_sd_10m, 1:6)

#Leaf Angle Mean
par(mfrow=c(2,3))
plot(Leaf_Angle_mean_10, 1:6)

gvlma(Leaf_Angle_mean_10)
## 
## Call:
## lm(formula = Leaf_Angle_mean ~ lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)  lulc_reduced149  lulc_reduced151  lulc_reduced155  
##          48.727            5.660            1.049            1.311  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Leaf_Angle_mean_10) 
## 
##                        Value p-value                Decision
## Global Stat        6.230e+00  0.1826 Assumptions acceptable.
## Skewness           2.686e+00  0.1012 Assumptions acceptable.
## Kurtosis           1.936e+00  0.1641 Assumptions acceptable.
## Link Function      5.268e-13  1.0000 Assumptions acceptable.
## Heteroscedasticity 1.608e+00  0.2048 Assumptions acceptable.
# Call:
# lm(formula = Leaf_Angle_mean ~ lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)  lulc_reduced149  lulc_reduced151  lulc_reduced155  
#          48.727            5.660            1.049            1.311  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Leaf_Angle_mean_10) 
# 
#                        Value p-value                Decision
# Global Stat        6.230e+00  0.1826 Assumptions acceptable.
# Skewness           2.686e+00  0.1012 Assumptions acceptable.
# Kurtosis           1.936e+00  0.1641 Assumptions acceptable.
# Link Function      5.268e-13  1.0000 Assumptions acceptable.
# Heteroscedasticity 1.608e+00  0.2048 Assumptions acceptable.

#Leaf Angle SD
par(mfrow=c(2,3))
plot(Leaf_Angle_sd_10, 1:6)

gvlma(Leaf_Angle_sd_10)
## 
## Call:
## lm(formula = Leaf_Angle_sd ~ tpi1000m + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)         tpi1000m  lulc_reduced149  lulc_reduced151  
##         24.7175          -4.1758          -6.5279          -0.9586  
## lulc_reduced155  
##         -3.1830  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Leaf_Angle_sd_10) 
## 
##                      Value p-value                Decision
## Global Stat        4.89145  0.2986 Assumptions acceptable.
## Skewness           0.48586  0.4858 Assumptions acceptable.
## Kurtosis           0.71337  0.3983 Assumptions acceptable.
## Link Function      3.67004  0.0554 Assumptions acceptable.
## Heteroscedasticity 0.02219  0.8816 Assumptions acceptable.
# Call:
# lm(formula = Leaf_Angle_sd ~ tpi1000m + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)         tpi1000m  lulc_reduced149  lulc_reduced151  lulc_reduced155  
#         24.7175          -4.1758          -6.5279          -0.9586          -3.1830  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Leaf_Angle_sd_10) 
# 
#                      Value p-value                Decision
# Global Stat        4.89145  0.2986 Assumptions acceptable.
# Skewness           0.48586  0.4858 Assumptions acceptable.
# Kurtosis           0.71337  0.3983 Assumptions acceptable.
# Link Function      3.67004  0.0554 Assumptions acceptable.
# Heteroscedasticity 0.02219  0.8816 Assumptions acceptable

#Diversity
par(mfrow=c(2,3))
plot(Diversity_10m, 1:6)

gvlma(Diversity_10m)
## 
## Call:
## lm(formula = div ~ DEM + northsouth + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)              DEM       northsouth  lulc_reduced149  
##         0.38439          0.09229         -0.04444          0.16987  
## lulc_reduced151  lulc_reduced155  
##         0.03923          0.15174  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Diversity_10m) 
## 
##                      Value  p-value                   Decision
## Global Stat        17.0809 0.001864 Assumptions NOT satisfied!
## Skewness            6.2070 0.012725 Assumptions NOT satisfied!
## Kurtosis            5.4838 0.019194 Assumptions NOT satisfied!
## Link Function       0.8535 0.355558    Assumptions acceptable.
## Heteroscedasticity  4.5366 0.033178 Assumptions NOT satisfied!
# Call:
# lm(formula = div ~ DEM + northsouth + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)              DEM       northsouth  lulc_reduced149  lulc_reduced151  
#         0.38439          0.09229         -0.04444          0.16987          0.03923  
# lulc_reduced155  
#         0.15174  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Diversity_10m) 
# 
#                      Value  p-value                   Decision
# Global Stat        17.0809 0.001864 Assumptions NOT satisfied!
# Skewness            6.2070 0.012725 Assumptions NOT satisfied!
# Kurtosis            5.4838 0.019194 Assumptions NOT satisfied!
# Link Function       0.8535 0.355558    Assumptions acceptable.
# Heteroscedasticity  4.5366 0.033178 Assumptions NOT satisfied!

#DBH Mean
par(mfrow=c(2,3))
plot(mvDBH_mean_10m, 1:6)

#DBH SD
par(mfrow=c(2,3))
plot(mvDBH_sd_10m, 1:6)

#Stand Density
par(mfrow=c(2,3))
plot(Stand_Density_10m, 1:6)

gvlma(Stand_Density_10m)
## 
## Call:
## lm(formula = Stand_Density ~ DEM + tpi100m + lulc_reduced, data = bighornData)
## 
## Coefficients:
##     (Intercept)              DEM          tpi100m  lulc_reduced149  
##          6.0015           1.2977           1.3722           3.0350  
## lulc_reduced151  lulc_reduced155  
##         -0.2019           1.5211  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Stand_Density_10m) 
## 
##                     Value p-value                Decision
## Global Stat        8.1492 0.08626 Assumptions acceptable.
## Skewness           2.6182 0.10564 Assumptions acceptable.
## Kurtosis           1.3827 0.23964 Assumptions acceptable.
## Link Function      3.3454 0.06739 Assumptions acceptable.
## Heteroscedasticity 0.8029 0.37023 Assumptions acceptable.
# Call:
# lm(formula = Stand_Density ~ DEM + tpi100m + lulc_reduced, data = bighornData)
# 
# Coefficients:
#     (Intercept)              DEM          tpi100m  lulc_reduced149  lulc_reduced151  lulc_reduced155 
#          6.0015           1.2977           1.3722           3.0350          -0.2019   1.5211  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Stand_Density_10m) 
# 
#                     Value p-value                Decision
# Global Stat        8.1492 0.08626 Assumptions acceptable.
# Skewness           2.6182 0.10564 Assumptions acceptable.
# Kurtosis           1.3827 0.23964 Assumptions acceptable.
# Link Function      3.3454 0.06739 Assumptions acceptable.
# Heteroscedasticity 0.8029 0.37023 Assumptions acceptable.


# Tree Height Mean
par(mfrow=c(2,3))
plot(Mean_Tree_Height10m, 1:6)

gvlma(Mean_Tree_Height10m)
## 
## Call:
## lm(formula = Mean_Tree_Height ~ northsouth + lulc_reduced, data = bighornDataTreeHeightsAndAge)
## 
## Coefficients:
##     (Intercept)       northsouth  lulc_reduced149  lulc_reduced151  
##          7.3446          -0.8025           5.3506           5.5459  
## lulc_reduced155  
##          6.9510  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Mean_Tree_Height10m) 
## 
##                       Value p-value                Decision
## Global Stat        1.052881  0.9017 Assumptions acceptable.
## Skewness           0.997379  0.3179 Assumptions acceptable.
## Kurtosis           0.036870  0.8477 Assumptions acceptable.
## Link Function      0.005801  0.9393 Assumptions acceptable.
## Heteroscedasticity 0.012832  0.9098 Assumptions acceptable.
# Call:
# lm(formula = Mean_Tree_Height ~ northsouth + lulc_reduced, data = bighornDataTreeHeightsAndAge)
# 
# Coefficients:
#     (Intercept)       northsouth  lulc_reduced149  lulc_reduced151  lulc_reduced155 
#          7.3446          -0.8025           5.3506           5.5459          6.9510   
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Mean_Tree_Height10m) 
# 
#                       Value p-value                Decision
# Global Stat        1.052881  0.9017 Assumptions acceptable.
# Skewness           0.997379  0.3179 Assumptions acceptable.
# Kurtosis           0.036870  0.8477 Assumptions acceptable.
# Link Function      0.005801  0.9393 Assumptions acceptable.
# Heteroscedasticity 0.012832  0.9098 Assumptions acceptable.

# Tree Height SD
par(mfrow=c(2,3))
plot(SD_Tree_Height10m, 1:6)

gvlma(SD_Tree_Height10m)
## 
## Call:
## lm(formula = Sd_Tree_Height ~ tpi100m + lulc_reduced, data = bighornDataTreeHeightsAndAge)
## 
## Coefficients:
##     (Intercept)          tpi100m  lulc_reduced149  lulc_reduced151  
##          5.0266          -0.8031          -1.0080           0.7790  
## lulc_reduced155  
##         -0.8551  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = SD_Tree_Height10m) 
## 
##                     Value p-value                   Decision
## Global Stat        4.9753 0.28984    Assumptions acceptable.
## Skewness           3.9299 0.04744 Assumptions NOT satisfied!
## Kurtosis           0.3955 0.52942    Assumptions acceptable.
## Link Function      0.3909 0.53182    Assumptions acceptable.
## Heteroscedasticity 0.2590 0.61080    Assumptions acceptable.
# Call:
# lm(formula = Sd_Tree_Height ~ tpi100m + lulc_reduced, data = bighornDataTreeHeightsAndAge)
# 
# Coefficients:
#     (Intercept)          tpi100m  lulc_reduced149  lulc_reduced151  lulc_reduced155
#          5.0266          -0.8031          -1.0080           0.7790    -0.8551  
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = SD_Tree_Height10m) 
# 
#                     Value p-value                   Decision
# Global Stat        4.9753 0.28984    Assumptions acceptable.
# Skewness           3.9299 0.04744 Assumptions NOT satisfied!
# Kurtosis           0.3955 0.52942    Assumptions acceptable.
# Link Function      0.3909 0.53182    Assumptions acceptable.
# Heteroscedasticity 0.2590 0.61080    Assumptions acceptable.

#Max Tree Height 
par(mfrow=c(2,3))
plot(Max_Tree_Height10m, 1:6)

gvlma(Max_Tree_Height10m)
## 
## Call:
## lm(formula = Max_Tree_Height ~ aspect + tpi1000m + lulc_reduced, 
##     data = bighornDataTreeHeightsAndAge)
## 
## Coefficients:
##     (Intercept)           aspect         tpi1000m  lulc_reduced149  
##         15.8762           0.8112          -1.1283           2.3294  
## lulc_reduced151  lulc_reduced155  
##          5.5133           6.5971  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Max_Tree_Height10m) 
## 
##                      Value p-value                Decision
## Global Stat        4.26114 0.37182 Assumptions acceptable.
## Skewness           0.46840 0.49372 Assumptions acceptable.
## Kurtosis           0.85119 0.35622 Assumptions acceptable.
## Link Function      2.91962 0.08751 Assumptions acceptable.
## Heteroscedasticity 0.02193 0.88228 Assumptions acceptable.
# Call:
# lm(formula = Max_Tree_Height ~ aspect + tpi1000m + lulc_reduced, 
#     data = bighornDataTreeHeightsAndAge)
# 
# Coefficients:
#     (Intercept)           aspect         tpi1000m  lulc_reduced149    lulc_reduced151  lulc_reduced155 
#         15.8762           0.8112          -1.1283           2.3294     5.5133           6.5971 
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Max_Tree_Height10m) 
# 
#                      Value p-value                Decision
# Global Stat        4.26114 0.37182 Assumptions acceptable.
# Skewness           0.46840 0.49372 Assumptions acceptable.
# Kurtosis           0.85119 0.35622 Assumptions acceptable.
# Link Function      2.91962 0.08751 Assumptions acceptable.
# Heteroscedasticity 0.02193 0.88228 Assumptions acceptable.

#Average Stand Age
par(mfrow=c(2,3))
plot(Average_Stand_Age10m, 1:6)

gvlma(Average_Stand_Age10m)
## 
## Call:
## lm(formula = Average_Stand_Age ~ aspect + tpi1000m + lulc_reduced, 
##     data = bighornDataTreeHeightsAndAge)
## 
## Coefficients:
##     (Intercept)           aspect         tpi1000m  lulc_reduced149  
##          74.491            5.211           -6.481          -19.671  
## lulc_reduced151  lulc_reduced155  
##           1.676          -18.959  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = Average_Stand_Age10m) 
## 
##                      Value  p-value                   Decision
## Global Stat        15.3846 0.003966 Assumptions NOT satisfied!
## Skewness            4.6874 0.030385 Assumptions NOT satisfied!
## Kurtosis            1.2152 0.270299    Assumptions acceptable.
## Link Function       0.8314 0.361869    Assumptions acceptable.
## Heteroscedasticity  8.6506 0.003270 Assumptions NOT satisfied!
# Call:
# lm(formula = Average_Stand_Age ~ aspect + tpi1000m + lulc_reduced, 
#     data = bighornDataTreeHeightsAndAge)
# 
# Coefficients:
#     (Intercept)           aspect         tpi1000m  lulc_reduced149  lulc_reduced151  lulc_reduced155 
#          74.491            5.211           -6.481          -19.671         l1.676          -18.959     
# 
# ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
# USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
# Level of Significance =  0.05 
# 
# Call:
#  gvlma(x = Average_Stand_Age10m) 
# 
#                      Value  p-value                   Decision
# Global Stat        15.3846 0.003966 Assumptions NOT satisfied!
# Skewness            4.6874 0.030385 Assumptions NOT satisfied!
# Kurtosis            1.2152 0.270299    Assumptions acceptable.
# Link Function       0.8314 0.361869    Assumptions acceptable.
# Heteroscedasticity  8.6506 0.003270 Assumptions NOT satisfied!

Gap Fraction Plots

#GapFraction_mean_10m<-lm(Gap_Fraction_mean~ TRI+ tpi1000m+lulc_reduced,data=bighornData)
GF.TRI<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=TRI, y=Gap_Fraction_mean, color=Gap_Fraction_mean)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Topographic Roughness Index", y ="Gap Fraction Mean",
       title ="Gap Fraction Mean vs Topographic Roughness IndexI", color ="Gap Fraction Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
GF.TRI

GF.100TPI<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi100m, y=Gap_Fraction_mean, color=Gap_Fraction_mean)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "100m TPI", y ="Gap Fraction Mean",
       title ="Gap Fraction Mean vs 100m TPI", color ="Gap Fraction Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
GF.100TPI

GF.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Gap_Fraction_mean, color=Gap_Fraction_mean)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Landcover Type", y ="Gap Fraction Mean",
       title ="Gap Fraction Mean vs Landcover Type", color ="Gap Fraction Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
GF.LULC

Gap_Fraction_mean_Car<-car::avPlots(model = GapFraction_mean_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=TRUE,
             main=("Gap Fraction Added-Variable Plots")) 

#GapFraction_sd_10m<-lm(Gap_Fraction_sd~  DEM, data=bighornData)
GFsd.DEM<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=Gap_Fraction_sd, color=Gap_Fraction_sd)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Elevation", y ="Gap Fraction Sd",
       title ="Gap Fraction Sd vs Elevation", color ="Gap Fraction Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
GFsd.DEM

Gap_Fraction_sd_Car<-car::avPlots(model = MVGapFraction_sd_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Gap Fraction Added-Variable Plot"))

PAR LAI Plots

#PAR LAI mean
#p-value: 0.1319
#PAR_LAI_mean_10m<-lm(PAR_LAI_mean~ DEM+tpi10m+lulc_reduced,data=bighornData)

PARLAI.DEM<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=PAR_LAI_mean, color=PAR_LAI_mean)) +
  scale_color_viridis_c(option="viridis")+
  labs(x = "Elevation", y ="PAR Leaf Area Index Mean",
       title ="PAR Leaf Area Index Mean vs Elevation", color ="PAR LAI Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
PARLAI.DEM

PARLAI.TPI10<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi10m, y=PAR_LAI_mean, color=PAR_LAI_mean)) +
  scale_color_viridis_c(option="viridis")+
  labs(x = "10m TPI", y ="PAR Leaf Area Index Mean",
       title ="PAR Leaf Area Index Mean vs 10m TPI", color ="PAR LAI Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
PARLAI.TPI10

PARLAI.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=PAR_LAI_mean, color=PAR_LAI_mean)) +
  scale_color_viridis_c(option="viridis")+
  labs(x = "Landcover Type", y ="PAR Leaf Area Index Mean",
       title ="PAR Leaf Area Index Mean vs Landcover Type", color ="PAR LAI Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
PARLAI.LULC

PAR_LAI_CAR<-car::avPlots(model = mvPAR_LAI_mean_10m, 
             col=carPalette()[7],
             col.lines=carPalette()[2], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=TRUE,
             main=("PAR LAI Added-Variable Plots")) 

#PAR LAI sd
#p-value: 0.02881
#PAR_LAI_sd_10m<-lm(PAR_LAI_sd~ eastwest, data=bighornData)
PARLAIsd.EW<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=eastwest, y=PAR_LAI_sd, color=PAR_LAI_sd)) +
  scale_color_viridis_c(option="viridis")+
  labs(x = "East-West Aspect", y ="PAR Leaf Area Index Sd",
       title ="PAR Leaf Area Index Sd vs East-West Aspect", color ="PAR LAI Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
PARLAIsd.EW

PAR_LAI.sd_CAR<-car::avPlots(model = PAR_LAI_sd_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("PAR LAI Added-Variable Plot"))

PAR Average Plots

#PAR Average mean
#p-value: 0.008254
#PAR_Average_mean_10m<-lm(PAR_Average_mean~  lulc_reduced,data=bighornData)

PARAvg.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=PAR_Average_mean, color=PAR_Average_mean)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Landcover Type", y ="PAR Average Mean",
       title ="PAR Average Mean vs Landcover Type", color ="PAR Average Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
PARAvg.LULC

PAR_Avg_CAR<-car::avPlots(model = PAR_Average_mean_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("PAR Average mean Added-Variable Plots")) 

#PAR Average SD
#p-value: 0.01623
#PAR_Average_sd10m<-lm(PAR_Average_sd~ TRI+slope+ northsouth,data=bighornData)
PARAvg.sd.TRI<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=TRI, y=PAR_Average_sd, color=PAR_Average_sd)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Topographic Roughness Index", y ="PAR Average Sd",
       title ="PAR Average Sd vs Topographic Roughness Index", color ="PAR Average Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
PARAvg.sd.TRI

PARAvg.sd.slope<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=slope, y=PAR_Average_sd, color=PAR_Average_sd)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Slope Anlge", y ="PAR Average Sd",
       title ="PAR Average Sd vs Slope Angle", color ="PAR Average Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
PARAvg.sd.slope

PARAvg.sd.NS<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=northsouth, y=PAR_Average_sd, color=PAR_Average_sd)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "North-South Aspect", y ="PAR Average Sd",
       title ="PAR Average Sd vs North-South Aspect", color ="PAR Average Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
PARAvg.sd.NS

PARAvg_Car<-car::avPlots(model = mvPAR_Average_sd10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("PAR Average sd Added-Variable Plots"))

Canopy Density Plots

#Canopy Density Mean 
#p-value: 8.287e-06
#Canopy_Density_mean_10m<-lm(Canopy_Density_mean~  DEM+ lulc_reduced,data=bighornData)

CD.DEM <-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=Canopy_Density_mean, color=Canopy_Density_mean)) +
  scale_color_viridis_c(option="rocket")+
  labs(x = "Elevation", y ="Canopy Density Mean",
       title ="Canopy Density Mean vs Elevation", color ="Canopy Density Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
CD.DEM 

CD.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Canopy_Density_mean, color=Canopy_Density_mean)) +
  scale_color_viridis_c(option="rocket")+
  labs(x = "Landcover Types", y ="Canopy Density Mean",
       title ="Canopy Density Mean vs Landcover Types", color ="Canopy Density Mean") + 
  theme(plot.title = element_text(hjust = 0.5))
CD.LULC

CD_CAR<-car::avPlots(model = Canopy_Density_mean_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Canopy Density mean Added-Variable Plots")) 

#Canopy Density SD
#5.175e-05
#Canopy_Density_sd_10m<-lm(Canopy_Density_sd~ DEM+tpi100m+ lulc_reduced, data=bighornData)
CDs.DEM <-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=Canopy_Density_sd, color=Canopy_Density_sd)) +
  scale_color_viridis_c(option="rocket")+
  labs(x = "Elevation", y ="Canopy Density Sd",
       title ="Canopy Density Sd vs Elevation", color ="Canopy Density Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
CDs.DEM 

CDs.TPI100<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi100m, y=Canopy_Density_sd, color=Canopy_Density_sd)) +
  scale_color_viridis_c(option="rocket")+
  labs(x = "100m TPI", y ="Canopy Density Sd",
       title ="Canopy Density Sd vs 100m TPI", color ="Canopy Density Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
CDs.TPI100

CDs.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Canopy_Density_sd, color=Canopy_Density_sd)) +
  scale_color_viridis_c(option="rocket")+
  labs(x = "Landcover Types", y ="Canopy Density Sd",
       title ="Canopy Density Sd vs Landcover Types", color ="Canopy Density Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
CDs.LULC

CDs_CAR<-car::avPlots(model = mvCanopy_Density_sd_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Canopy Density sd Added-Variable Plots"))

Leaf Angle Plots

#Leaf Angle mean
#p-value: 0.5935
#Leaf_Angle_mean_10<-lm(Leaf_Angle_mean~ lulc_reduced, data=bighornData)

LA.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Leaf_Angle_mean, color=Leaf_Angle_mean)) +
  scale_color_viridis_c(option="plasma")+
  labs(x = "Landcover Types", y ="Leaf Angle Mean",
       title ="Leaf Angle Mean vs Landcover Types", color ="Leaf Angle Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
LA.LULC

LA_CAR<-car::avPlots(model = Leaf_Angle_mean_10, 
             col=carPalette()[1],
             col.lines=carPalette()[1], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Leaf Angle mean Added-Variable Plots")) 

#Leaf Angle sd
#p-value: 0.001394
#Leaf_Angle_sd_10<-lm(Leaf_Angle_sd~ tpi1000m+lulc_reduced,data=bighornData)
LAs.TPI1000<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi1000m, y=Leaf_Angle_sd, color=Leaf_Angle_sd)) +
  scale_color_viridis_c(option="plasma")+
  labs(x = "1000m TPI", y ="Leaf Angle Sd",
       title ="Leaf Angle Sd vs 1000m TPI", color ="Leaf Angle Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
LAs.TPI1000

LAs.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Leaf_Angle_sd, color=Leaf_Angle_sd)) +
  scale_color_viridis_c(option="plasma")+
  labs(x = "Land Cover Types", y ="Leaf Angle Sd",
       title ="Leaf Angle Sd vs Land Cover Types", color ="Leaf Angle Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
LAs.LULC

LA.sd_CAR<-car::avPlots(model = Leaf_Angle_sd_10, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Leaf Angle sd Added-Variable Plots"))

Diversity Plots

#Diversity
# p-value: 0.0009821
#Diversity_10m<-lm(div~DEM+northsouth+lulc_reduced, data=bighornData)

div.DEM<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=div, color=div)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Elevation", y ="Diversity",
       title ="Diversity vs Elevation", color ="Diversity") + 
  theme(plot.title = element_text(hjust = 0.5))
div.DEM

div.NS<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=northsouth, y=div, color=div)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "North-South Aspect", y ="Diversity",
       title ="Diversity vs North-South Aspect", color ="Diversity") + 
  theme(plot.title = element_text(hjust = 0.5))
div.NS

div.lulc<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=div, color=div)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Landcover Types", y ="Diversity",
       title ="Diversity vs Landcover Types", color ="Diversity") + 
  theme(plot.title = element_text(hjust = 0.5))
div.lulc

div_Car<-car::avPlots(model = Diversity_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Diversity Added-Variable Plots"))

DBH Plots

# Average DBH
#p-value: 0.0009843
#DBH_mean_10m<-lm(Average_DBH~ tpi100m+ lulc_reduced,data=bighornData)
DBH_Avg.TPI100<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi100m, y=Average_DBH, color=Average_DBH)) +
  scale_color_viridis_c(option="cividis")+
  labs(x = "100m TPI", y ="Average DBH",
       title ="Average Diameter Breast Height vs 100m TPI", color ="Average DBH") + 
  theme(plot.title = element_text(hjust = 0.5))
DBH_Avg.TPI100

DBH_Avg.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Average_DBH, color=Average_DBH)) +
  scale_color_viridis_c(option="cividis")+
  labs(x = "Landcover Types", y ="Average DBH",
       title ="Average Diameter Breast Height vs Landcover Types", color ="Average DBH") + 
  theme(plot.title = element_text(hjust = 0.5))
DBH_Avg.LULC

DBH_Avg_CAR<-car::avPlots(model = mvDBH_mean_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("DBH mean Added-Variable Plots")) 

#SD DBH
#p-value: 0.001216
#DBH_sd_10m<-lm(DBH_Sd~ tpi1000m+ lulc_reduced, data=bighornData)
dbh.sd.tpi1000<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi1000m, y=DBH_Sd, color=DBH_Sd)) +
  scale_color_viridis_c(option="cividis")+
  labs(x = "1000m TPI", y ="DBH Sd",
       title ="Sd of Diameter Breast Height vs 1000m TPI", color ="DBH Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
dbh.sd.tpi1000

dbh.sd.LULC<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=DBH_Sd, color=DBH_Sd)) +
  scale_color_viridis_c(option="cividis")+
  labs(x = "Land Cover Types", y ="DBH Sd",
       title ="Sd of Diameter Breast Height vs Land Cover Types", color ="DBH Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
dbh.sd.LULC

lmDBH_sd_10m<-glm(DBH_Sd~ tpi1000m+ lulc_reduced, 
                family=gaussian(),
                data=bighornData) 

dbh.sd.CAR<-car::avPlots(model = lmDBH_sd_10m,
             col=carPalette()[3],
             col.lines=carPalette()[5],
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("DBH sd Added-Variable Plots"))

Stand Density Plots

#Stand Density
#p-value: 1.806e-06
#Stand_Density_10m<-lm(Stand_Density~ DEM+ tpi100m+ lulc_reduced, data=bighornData)

den.dem<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=DEM, y=Stand_Density, color=Stand_Density)) +
  scale_color_viridis_c(option="mako")+
  labs(x = "Elevation", y ="Stand Density",
       title ="Stand Density vs Elevation", color ="Stand Density") + 
  theme(plot.title = element_text(hjust = 0.5))
den.dem

den.tpi100<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=tpi100m, y=Stand_Density, color=Stand_Density)) +
  scale_color_viridis_c(option="mako")+
  labs(x = "100m TPI", y ="Stand Density",
       title ="Average Tree Height vs 100m TPI", color ="Stand Density") + 
  theme(plot.title = element_text(hjust = 0.5))
den.tpi100

den.lulc<-ggplot(bighornDataNS)+
  geom_point(mapping=aes(x=lulc_reduced, y=Stand_Density, color=Stand_Density)) +
  scale_color_viridis_c(option="mako")+
  labs(x = "Landcover Types", y ="Stand Density",
       title ="Stand Density vs Landcover Types", color ="Stand Density") + 
  theme(plot.title = element_text(hjust = 0.5))
den.lulc

den.car<-car::avPlots(model = Stand_Density_10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Stand Density Added-Variable Plots"))

Tree Height Plots

#Mean Tree Height
#p-value: 0.0199
#Mean_Tree_Height10m<-lm(Mean_Tree_Height~northsouth+lulc_reduced, data=bighornDataTreeHeightsAndAge)

avgTH.NS<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=northsouth, y=Mean_Tree_Height, color=Mean_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "North-South Aspect", y ="Average Tree Height",
       title ="Average Tree Height vs North-South Aspect", color ="Average Tree Height") + 
  theme(plot.title = element_text(hjust = 0.5))
avgTH.NS

avgTH.LULC<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=lulc_reduced, y=Mean_Tree_Height, color=Mean_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Land Cover Types", y ="Average Tree Height",
       title ="Average Tree Height vs Land Cover Types", color ="Average Tree Height") + 
  theme(plot.title = element_text(hjust = 0.5))
avgTH.LULC

avgTH.CAR<-car::avPlots(model = Mean_Tree_Height10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Average Tree Height Added-Variable Plots")) 

#SD Tree Height
#p-value: 0.006275
#SD_Tree_Height10m<-lm(Sd_Tree_Height~ tpi100m+lulc_reduced,  data=bighornDataTreeHeightsAndAge)

sd.th.tpi100<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=tpi100m, y=Sd_Tree_Height, color=Sd_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "100m TPI", y ="Tree Height Sd",
       title ="Tree Height Sd vs 100m TPI", color ="Tree Height Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
sd.th.tpi100

sd.th.lulc<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=lulc_reduced, y=Sd_Tree_Height, color=Sd_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Land Cover Types", y ="Tree Height Sd",
       title ="Tree Height Sd vs Aspect", color ="Tree Height Sd") + 
  theme(plot.title = element_text(hjust = 0.5))
sd.th.lulc

sd.thCAR<-car::avPlots(model = SD_Tree_Height10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Tree Height sd Added-Variable Plots")) 

#Max Tree Height (northsouth aspect)
#p-value: 0.00969
#Max_Tree_Height10m<-lm(Max_Tree_Height~ aspect+ tpi1000m+lulc_reduced,data=bighornDataTreeHeightsAndAge)
max.th.aspect<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=aspect, y=Max_Tree_Height, color=Max_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Aspect", y ="Max Tree Height",
       title ="Maximum Tree Heights vs Aspect", color ="Maximum Tree Height") + 
  theme(plot.title = element_text(hjust = 0.5))
max.th.aspect

max.th.tpi1000<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=tpi1000m, y=Max_Tree_Height, color=Max_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "1000m TPI", y ="Max Tree Height",
       title ="Maximum Tree Heights vs 1000m TPI", color ="Maximum Tree Height") + 
  theme(plot.title = element_text(hjust = 0.5))
max.th.tpi1000

max.th.lulc<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=lulc_reduced, y=Max_Tree_Height, color=Max_Tree_Height)) +
  scale_color_viridis_c(option="magma")+
  labs(x = "Landcover Types", y ="Max Tree Height",
       title ="Maximum Tree Heights vs Landcover Type", color ="Maximum Tree Height") + 
  theme(plot.title = element_text(hjust = 0.5))
max.th.lulc

max.th.CAR<-car::avPlots(model = Max_Tree_Height10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Maximum Tree Height Added-Variable Plots"))

Stand Age Plots

#Average Stand Age 
#p-value: 0.0001305
#Average_Stand_Age10m<-lm(Average_Stand_Age~aspect+ tpi1000m +lulc_reduced,  data=bighornDataTreeHeightsAndAge)
avg.sa.aspect<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=aspect, y=Average_Stand_Age, color=Average_Stand_Age)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "aspect", y ="Average_Stand_Age", title ="Average Stand Age vs Aspect", color ="Average_Stand_Age") +
  theme(plot.title = element_text(hjust = 0.5))
avg.sa.aspect

avg.sa.tpi1000<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=tpi1000m, y=Average_Stand_Age, color=Average_Stand_Age)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "1000m TPI", y ="Average Stand Age", title ="Average Stand Age vs 1000m TPI", color ="Average Stand Age") +
  theme(plot.title = element_text(hjust = 0.5))
avg.sa.tpi1000

avg.sa.lulc<-ggplot(NSbighornDataTreeHeightsAndAge)+
  geom_point(mapping=aes(x=lulc_reduced, y=Average_Stand_Age, color=Average_Stand_Age)) +
  scale_color_viridis_c(option="turbo")+
  labs(x = "Landcover Types", y ="Average Stand Age",
       title ="Average Stand Age vs Landcover Type", color ="Average Stand Age") + 
  theme(plot.title = element_text(hjust = 0.5))
avg.sa.lulc

avg.sa.CAR<-car::avPlots(model = Average_Stand_Age10m, 
             col=carPalette()[3],
             col.lines=carPalette()[5], 
             grid=TRUE,
             ellipse=FALSE,
             marginal.scale=FALSE,
             main=("Estimated Stand Age Added-Variable Plots"))